// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtDualMeshFilter.C
// ************************************************************************* //

#include <avtDualMeshFilter.h>

#include <vtkDataSet.h>
#include <vtkDataArray.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <vtkRectilinearGrid.h>

// ****************************************************************************
//  Method: avtDualMeshFilter constructor
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed May 7 15:59:34 PST 2008
//
// ****************************************************************************

avtDualMeshFilter::avtDualMeshFilter()
{
    actualVar  = "";
    actualMode = DualMeshAttributes::Auto;
}


// ****************************************************************************
//  Method: avtDualMeshFilter destructor
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed May 7 15:59:34 PST 2008
//
//  Modifications:
//
// ****************************************************************************

avtDualMeshFilter::~avtDualMeshFilter()
{
}


// ****************************************************************************
//  Method:  avtDualMeshFilter::Create
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed May 7 15:59:34 PST 2008
//
// ****************************************************************************

avtFilter *
avtDualMeshFilter::Create()
{
    return new avtDualMeshFilter();
}


// ****************************************************************************
//  Method:      avtDualMeshFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed May 7 15:59:34 PST 2008
//
// ****************************************************************************

void
avtDualMeshFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const DualMeshAttributes*)a;
}


// ****************************************************************************
//  Method: avtDualMeshFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtDualMeshFilter with the given
//      parameters would result in an equivalent avtDualMeshFilter.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed May 7 15:59:34 PST 2008
//
// ****************************************************************************

bool
avtDualMeshFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(DualMeshAttributes*)a);
}

// ****************************************************************************
//  Method: avtDualMeshFilter::UpdateDataObjectInfo
//
//  Purpose:
//      Indicates the zones no longer correspond to the original problem.
//
//  Programmer: Cyrus Harrison
//  Creation:   May 8, 2008
//
//  Modifications:
//    Cyrus Harrison, Tue May 13 13:42:13 PDT 2008
//    Fixed mode setting for empy meshes. 
//
//    Brad Whitlock, Mon Apr  7 15:55:02 PDT 2014
//    Add filter metadata used in export.
//    Work partially supported by DOE Grant SC0007548.
//
// ****************************************************************************

void
avtDualMeshFilter::UpdateDataObjectInfo(void)
{
    actualMode = atts.GetMode();
    
    GetOutput()->GetInfo().GetValidity().InvalidateZones();
    GetOutput()->GetInfo().GetValidity().InvalidateZones();
    
    avtDataAttributes &in_atts = GetInput()->GetInfo().GetAttributes();
    if (!in_atts.ValidActiveVariable())
    {
        // We don't have enough information to resolve Auto mode.
        return;
    }
    
        
    if(atts.GetMode() == DualMeshAttributes::Auto)
    {
        actualMode = DualMeshAttributes::NodesToZones;
        // if the active variable is not a mesh, use it to determine
        // conversion mode
        avtDataAttributes &in_atts = GetInput()->GetInfo().GetAttributes();

        if (in_atts.ValidVariable(actualVar.c_str()) &&
            in_atts.GetCentering(actualVar.c_str()) == AVT_ZONECENT)
        {
            actualMode = DualMeshAttributes::ZonesToNodes;
        }
    }

    
    if( actualMode == DualMeshAttributes::NodesToZones)
        GetOutput()->GetInfo().GetAttributes().SetCentering(AVT_ZONECENT);
    else if ( actualMode == DualMeshAttributes::ZonesToNodes)
        GetOutput()->GetInfo().GetAttributes().SetCentering(AVT_NODECENT);

    GetOutput()->GetInfo().GetAttributes().AddFilterMetaData("DualMesh");
}

// ****************************************************************************
//  Method: avtDualMeshFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the DualMesh filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed May 7 15:59:34 PST 2008
//
//  Modifications:
//    Cyrus Harrison, Fri Feb 13 15:03:08 PST 2009
//    Fixed memory leak.
//
//    Hank Childs, Tue Feb  8 15:28:14 PST 2011
//    Add better support for 2D data ... z-coord wasn't being set.
//
//    Eric Brugger, Wed Jul 23 12:09:19 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

avtDataRepresentation *
avtDualMeshFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();

    // require Rectilinear Grid 
    int  mesh_type = in_ds->GetDataObjectType();
    if (mesh_type != VTK_RECTILINEAR_GRID)
    {
        EXCEPTION1(ImproperUseException,
                   "The Dual Mesh operator only works on rectilinear grids.");
    }
    
    vtkRectilinearGrid *rgrid  = (vtkRectilinearGrid *) in_ds;
    vtkRectilinearGrid *result = vtkRectilinearGrid::New();
    
    // shallow copy the field data
    result->GetFieldData()->ShallowCopy(rgrid->GetFieldData());
    
    int dims[3]  = {0,0,0};
    int rdims[3] = {1,1,1};
    
    rgrid->GetDimensions(dims);

    // see if we are expanding or contracting
    bool expand = actualMode != DualMeshAttributes::ZonesToNodes;

    vtkDataArray *x_coords = rgrid->GetXCoordinates();
    vtkDataArray *y_coords = rgrid->GetYCoordinates();
    vtkDataArray *z_coords = rgrid->GetZCoordinates();
            
    vtkDataArray *res_x_coords = NULL;
    vtkDataArray *res_y_coords = NULL;
    vtkDataArray *res_z_coords = NULL;
        
    
    if(expand)
        res_x_coords  = ExpandDual(x_coords);
    else        
        res_x_coords  = ContractDual(x_coords);
        
    rdims[0] = res_x_coords->GetNumberOfTuples();
        
    if( dims[1] > 1 )  // 2d & 3d
    {
        if(expand)
            res_y_coords  = ExpandDual(y_coords);
        else
            res_y_coords  = ContractDual(y_coords);

        rdims[1] = res_y_coords->GetNumberOfTuples();
    }
        
    if( dims[2] > 1 )  // 3d
    {
        if(expand)
            res_z_coords  = ExpandDual(z_coords);
        else
            res_z_coords  = ContractDual(z_coords);

        rdims[2] = res_z_coords->GetNumberOfTuples();
    }
        
    result->SetDimensions(rdims);
    
    result->SetXCoordinates(res_x_coords);

    if(dims[1]>1)
        result->SetYCoordinates(res_y_coords);
    else
        result->SetYCoordinates(rgrid->GetYCoordinates());
            
    if(dims[2]>1)
        result->SetZCoordinates(res_z_coords);
    else
        result->SetZCoordinates(rgrid->GetZCoordinates());
            
    if(res_x_coords)
        res_x_coords->Delete();
    if(res_y_coords)
        res_y_coords->Delete();
    if(res_z_coords)
        res_z_coords->Delete();

    if(expand)
    {   
        // change point data to cell data 
        result->GetCellData()->ShallowCopy(rgrid->GetPointData());
    }
    else
    {
        // change cell data to point data
        result->GetPointData()->ShallowCopy(rgrid->GetCellData());
    }
    
    avtDataRepresentation *out_dr = new avtDataRepresentation(result,
        in_dr->GetDomain(), in_dr->GetLabel());

    result->Delete();
    
    return out_dr;
}

// ****************************************************************************
//  Method: avtDualMeshFilter::ExpandDual
//
//  Purpose:
//      Calculate the new coordiantes in one dimension for zonal to nodal 
//      dual mesh conversion.
//
//  Arguments:
//      coords    Input coordinate array.
//
//  Returns:      The new coordiante array.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed May 7 15:59:34 PST 2008
//
//  Modifications:
//    Kathleen Biagas, Tue Oct 16 14:01:34 MST 2012
//    Preserve coordinate type.
//
// ****************************************************************************

vtkDataArray *
avtDualMeshFilter::ExpandDual(vtkDataArray *coords)
{
    // compute expanded dims (cell to point)
    double diff    = 0.0;
    double curr    = 0.0;
    int   npoints = coords->GetNumberOfTuples();
        
    vtkDataArray *res_coords = coords->NewInstance();
        
    for(int i=0; i< npoints; i++)
    {
        curr = coords->GetTuple1(i);
        if(i < npoints - 1)
            diff = (coords->GetTuple1(i+1) + curr) / 2.0 - curr;
        
        if(i == 0)
        {
            // do first point, use coord[0] - (coord[1] - coord[0])/2
            res_coords->InsertNextTuple1(curr - diff);    
        }
        res_coords->InsertNextTuple1(curr + diff);
    }

    return res_coords;
}

// ****************************************************************************
//  Method: avtDualMeshFilter::ContractDual
//
//  Purpose:
//      Calculate the new coordiantes in one dimension for nodal to zonal
//      dual mesh conversion.
//
//  Arguments:
//      coords    Input coordinate array.
//
//  Returns:      The new coordiante array.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed May 7 15:59:34 PST 2008
//
//    Kathleen Biagas, Thu Oct 18 10:01:46 PDT 2012
//    Preserve coordinate data type.
//
// ****************************************************************************

vtkDataArray *
avtDualMeshFilter::ContractDual(vtkDataArray *coords)
{
     // compute contracted dims ( point to cell ) 
    double curr    = 0.0;
    int    npoints = coords->GetNumberOfTuples();
        
    vtkDataArray *res_coords = coords->NewInstance();
        
    for(int i=0;i< npoints - 1;i++)
    {
        curr = (coords->GetTuple1(i+1) + coords->GetTuple1(i)) / 2.0;
        res_coords->InsertNextTuple1(curr);
    }
    return res_coords;
}

// ****************************************************************************
//  Method: avtDualMeshFilter::ExamineContract
//
//  Purpose:
//      Gets the name of the active variable, which is needed if the "Auto"
//      conversion mode is enabled. 
//
//  Programmer: Cyrus Harrison
//  Creation:   May 8, 2008
//
// ****************************************************************************

void 
avtDualMeshFilter::ExamineContract(avtContract_p contract)
{
    avtDatasetToDatasetFilter::ExamineContract(contract);
    actualVar = contract->GetDataRequest()->GetVariable();
}


// ****************************************************************************
// Method: avtDualMeshFilter::ModifyContract
//
// Purpose: 
//   Modify the data request in the contract.
//
// Arguments:
//   spec : The pipeline specification.
//
// Returns:    A new pipeline specification.
//
// Programmer: Brad Whitlock
// Creation:   Wed Aug 15 12:23:39 PDT 2012
//
// Modifications:
//
// ****************************************************************************

avtContract_p
avtDualMeshFilter::ModifyContract(avtContract_p spec)
{
    avtContract_p rv = new avtContract(spec);

    // Indicate that we want to ignore missing data instead of removing it.
    rv->GetDataRequest()->IgnoreMissingData();

    return rv;
}
